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Abstract. Using the shift-operator technique, a compact formula for the Fourier 
transform of a product of two Slater-type orbitals located on different atomic centers 
is derived. The result is valid for arbitrary quantum numbers and was found to be 
numerically stable for a wide range of geometrical parameters and momenta. Details 
of the implementation are presented together with benchmark data for representative 
integrals. We also discuss the assets and drawbacks of alternative algorithms available 
and analyze the numerical efficiency of the new scheme. 
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1. Introduction 

The electron-electron interaction as quantified in Coulomb or exchange integrals is at the 
heart of every quantum mechanical treatment of condensed matter. Due to the simple 
structure of the Coulomb operator in reciprocal space, Fourier transform techniques 
allow for the transformation of the double integral over real space into a compact single 
integral in momentum space. Let us consider a typical two-electron repulsion integral 
as an example 



If we denote the Fourier transform of the product of orbitals M (r — R^) and </v(r — Rg) 
by 4>^ u (Ra, Rb, q), we may very schematically write (thorough definitions follow later) 



Methods along these lines date back at least to Bonham et al. pfl [2] and are the de-facto 
standard for systems with translational symmetry. For these boundary conditions, plane 
waves are the most natural type of basis functions, thanks also to their trivial behavior 
under Fourier transformation. 

In recent years, correlated electronic structure methods like the GW approximation 
of Hedin [3J, which were originally developed in the context of band structure 
calculations, are now also applied to systems where translational symmetry is broken. 
Examples of this kind include super lattices, defects, surfaces and even atoms or 
molecules jU EJ EJ [7J [H] . Obviously, atomic orbital basis sets are more appropriate in these 
situations and reduce the number of basis functions needed to achieve a certain accuracy 
considerably. Among the localized basis sets, Slater and Gaussian-type orbitals are the 
most prominent. The latter have the advantage that the Fourier transform of orbital 
products is relatively easy to obtain, while usually much less Slater than Gaussian-type 
orbitals are required to represent atomic or molecular electron densities. In the context 
of the Fourier transform methods mentioned above, the choice of Slater versus Gaussian 
is hence intimately connected with the ability to perform the Fourier transform of basis 
function products efficiently 

A direct numerical quadrature of the three dimensional integral over reciprocal 
space using Fast Fourier Transform techniques is not advisable due to high memory 
consumption, low computational speed and very limited numerical accuracy. As an 
alternative, Slater functions may be fitted to a fixed linear combination of Gaussian-type 
orbitals like done for example in the popular Pople basis sets often employed in quantum 
chemistry P, [TO] - However, in this case the electronic structure calculation could 
have been performed entirely in terms of Gaussians in the first place with additional 
variational freedom. Several attempts to directly perform the intricate integration of 
Slater products analytically are documented in the literature. While the mentioned 
original work of Bonham [H [2] was restricted to s-type functions, Bentley and Stewart 
[TT] derived an expression for arbitrary angular momentum states involving an infinite 
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series. Later, Junker [T2] obtained a result in terms of finite sums and one-dimensional 
numerical integrations, while Straton [T3J, [T3] was able to provide general formulas for 
the Fourier transform of the product of more than two orbitals. 

This earlier work meets some but not all desired properties of a general solution. 
The latter include the validity for arbitrary quantum numbers, the possibility for a 
straightforward implementation on a computer as well as high numerical efficiency and 
stability. Moreover, the solution should be amenable to partial wave analysis, in order 
to allow for an efficient evaluation of two-electron integrals in a second step. This point 
is maybe not so important for periodic systems, where quadrature of the remaining 
integral over reciprocal space in (T5]) may be accomplished by summation over special 
k-point meshes [15J , involving only a small number of integrand evaluations. For finite 
systems however, this point becomes crucial. Here it should be mentioned that the 
efficient evaluation of Fourier transforms is only a first step in a fast computation of 
multi-centre integrals. 

An approach which combines most of the above mentioned merits was proposed by 
Trivedi and Steinborn [16] . The authors provide formulas for the Fourier transform over 
products of so-called B-functions. These B-functions can be transformed into Slater-type 
orbitals without loss of generality and accuracy. Subsequently, this technique was used 
by Grotendorst and Steinborn [17] to evaluate a variety of multi-centre integrals required 
in electronic structure calculations. Alternative representations of the transforms are 
given in [18] and in the dissertation of Homeier [19], which also contains a deeper 
discussion of the B-function formalism together with numerical results and benchmark 
data. 

In this work, we propose an alternative to the Trivedi-Steinborn formula, which is 
directly formulated in terms of Slater-type orbitals. The derivation is based on the shift- 
operator technique, which is discussed in the next section. The approach may be seen 
as a generalization of a recent result for overlap integrals [20] and meets all important 
criteria established above. 

The more general aim of this contribution is to facilitate the utilization of Slater- 
type orbitals in the simulation of periodic and quasi-periodic systems. Currently, only a 
very limited number of codes employs this kind of basis set [21], [22] , due to the apparent 
difficulties in the numerical implementation. The development of adapted algorithms 
is therefore of key importance in order to unveil the well known benefits of Slater-type 
orbitals in atomistic calculations. 

2. Definitions and the shift-operator approach 

We consider real unnormalized Slater-type orbitals (STO) of the form: 




(3) 
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where n = n — I in terms of the principle quantum number n. The regular harmonics 
zf 1 are related to the more familiar real spherical harmonics Yj m : 

z?{t) = r l Y lm (v) (4) 

Y lm (r) = (-l) OT P™(cos0) cos(m0) m > 

yj_ ro (f) = (-l) m P, ro (cos0) sin(m0) m < 0, 
where P z m (cos#) denote associated Legendre polynomials as defined in [23]. Square 
normalized STO Xmm are readily obtained as: 



Xnlm 



(2Z + 1) (/ - |m|)! (2C) Mi+1 

Xnim- (5) 



\J 27r(l + 5 m0 ) (i + H)! (2n + 2Z)! 
The idea of the shift-operator approach is to evaluate the desired integral of interest, 
e.g., overlap or two-electron repulsion, first for the simplest STO of s-symmetry, for 
which an quadrature is often relatively easy. In a second step, the quantum numbers 
are then raised by operators that involve derivatives with respect to parameters of the 
integral, like decay constants and inter-center distance. The benefit of such raising and 
lowering operators in the solution of molecular integrals was recognized quite early and 
exploited by various authors [211 [251 126| [271 128| [29] . 

In this approach, a STO centered at R/ as a function of 17 = r — Rf may be written 

as: 

Xnlm ( ri ,() = nuvi)^-, (6) 

with V/ denoting the vector (d/dX I ,d/dY I ,d/dZ I ) and 

^(V,) = ^(V f )(-|)"(-i|)'. (7) 

A detailed discussion of the properties of ^ m (V/) and related differential operators is 
provided in a recent review by Weniger [30] . 

The form in (jTJ) is now used to construct the Fourier transform of two-center STO 
products which are in the focus of this work: 

OS(k,Ci,C2,R/,Rj)= / dre^ X n lh m^iXi)Xn 2 i 2m2 (rjX2) (8) 

= ^(V/^UVj) / C ?re to X ooo(rj,Ci)Xooo(r J ,C2). (9) 

The shift-operator approach is applicable if the basic integrals (Iqqq in our case) have 
a closed form which can be easily differentiated with respect to the outer parameters. 
The next section shows that this is indeed the case for the present Fourier transform. 

3. The basic integral 

As shown by Rico and co-workers [31], the product of two s-type STO can be expressed 
as an one-dimensional integral which is suitable for further manipulations: 



1 f 1 -2 > 

Xooo(r/, Ci)Xooo(r j, Ca) = - / du M 1 ~ u )} 2 C« k -i I C« 
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with R = Rj — R/, H u = uHj + (1 — it)Rj, r u = r — R u , = (fu + £f(l — u), and 
fc„(a;) = x v K v [x\ where K v (x) is the modified Bessel function of the second kind, often 
also called McDonald function. 

After insertion of (jTOl) into jSJ) and change of the integration variable to r n , the 
angular integration is readily performed by expanding the exponential in partial waves 



e4kr = sj^ r jy^^m^os9)J l+l {kr), = Z(k,r), (11) 

and using the fact that the remainder of the integrand has s-symmetry. The result 
reads: 

Ok,Ci,C2,R/,Rj) = r C du^[u(l-u)]-hl (12) 



o 



k 

x / dr u r u sm(kr u )k-i I ( u < 



o 



W + ^L^ . (13) 
\ u(l — n) 

The remaining radial integral is known [32] . which leads to the final result for the basic 
integral: 

/ ° °(k,Ci,C2,R/,Rj) = V^Rj^due 1 ™^ (Ry/Q + k*u(l-u)) ■ (14) 

For vanishing momentum transfer this formula reduces to the known result for the 
corresponding overlap integral as given, e.g., by Ema et al. [20J. For this special 
case, the pending integral may be solved analytically and is related to confluent 
hypergeometric functions. In general however, an evaluation based on numerical 
integration is unavoidable at this point. 

4. Transforms for higher quantum numbers 

Using the shift-operator approach, Fourier transforms for higher quantum numbers may 
now be written as: 

C£ 2 1 (k,Ci,C 2 ,R/,Rj) = z^iV^iVj) £ due ik ^Kl l l 1 2 (k,Ci,C2,R,u) (15) 



with 



d V 1 ( id V 1 



dCiJ V CidCi 

d V 2 ( id ' 



ki[RjQ + k 2 u(l-u)), (16) 



where, as shown in the Appendix A the derivatives with respect to the decay constants 
are relatively easy to perform. The action of the solid harmonics on the integral is more 
involved and requires further consideration. We proceed by introducing the equality 

^(V)(/s) = E E £ df m , mll {ztAV) f) {zf {V) g) (17) 

l'=0m'=-(l-l') m"=-l> 
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for arbitrary functions /(R),g(R). This relation is proven in the Appendix B using 



the Leibniz rule for the differentiation of products together with the completeness and 
orthogonality relations of spherical harmonics. Alternative proofs were given by Dunlap 
[33] and Weniger [30] . 

Values for the coefficients cZ^, m « can be obtained by straightforward differentiation 
for small quantum numbers. In general, the use of symbolic computation software 
allows the determination once and for all. Special cases include dQ™, m „ = 6 m / m 5 m 'i and 

^Im'm" = $m'0$m"m- 



Applying (fT7|) to (1T51) . we arrive at 

h h—l'x l'i 

C£( k 'CliC2>R/>Rj) = E ih ~ l1 E E ^llml'mx" 

l' 1= mi'=-(2i-ii) mi"=-/J 

l 2 l2—l' 2 l'i -1 

E ^ E 



X 



l' 2 =0 



U l' 2 m 2 'm 2 " I 

-(l 2 -l' 2 )m 2 "=-l' 2 



du 



u 



h-i'm 



u 



h-l' 2 jmi' 



zT^K-i^v^ 



X 



u 



[zf\V I )zf\V J )h^ 2 (k,CiX2,R,' 
where we used the homogeneity of regular harmonics and the fact that plane waves are 
eigenfunctions of the momentum operator. The remaining derivation parallels the work 
of Ema et al. [20] on overlap integrals and we will follow the nomenclature used there as 



close as possible to facilitate comparison. Since Zi™^* ^ n ^ ne ^ as ^ ^ me °f depends only 
on the norm of R, the following theorem may be applied which goes back to Hobson 

^ f(R) (19) 



C(V/)C(v,)/(^) 



L< <^—k 



-i) h Y — 



k=0 



k=0 



RdR / 

Him) W 



Here L< = min(Zi, l 2 ) and the V l k imihm2 are given by 



V: 



l\m\l 2 m 2 



)fc L< 



Z!r(Zi + h - 1 + 3/2)R 



2l-2k 



hm 1 l 2 m 2 m /p\ 
•w. I -TV- J 



+l 2 -2lm^h+l 2 -2l\ 



(20) 



(21) 



, _ fc (Z - fc)!r(Zx + Z 2 - Z - + 3/2) ,„ 

where the coefficients cjjT/ Z -Tzm are directly related to real Gaunt coefficients (For a 
detailed derivation of (fl9j) to ( J2TI) see Appendix C ). 

Next we define the quantity S$f 2h (this is a generalization of <5™ liin2 ' 2 i n the work 
of Ema et al. [20]), which is further discussed in the Appendix A 
~^W^ k;Ci)C2)R/)Rj) 



due ik ^u l2 - l Hl 



RdR 

n 2 



K±(kXi,(2,R,u) 



X 



i=L ry+l Ji=L ^j 
v 7 ^^ 1 cZwe 4kR "^(l - u)"jfe a + k 2 u(l - 1 



(22) 
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with L^J denoting the integer part of r and 

H = li + l 2 -l' 2 + i; v = l 1 + l 2 -l[+j\ a = -l/2-l 1 -l 2 + l' 1 + l' 2 -i-j-k 

(2()n(2i - n)\(n - 1 ' 

With these definitions we reach the main result of this work: 

CS 2 1 ( k »Ci,c 2 ,R/,R,) = E(-i)^ E €-i[( k ) E rf & mi »x (24) 

Z^=0 mi'=-(«i-Z;) mi"=-/i 

E E 2 C4( k ) E ^ w E ^""^n^^j/r^k, Ci, C2, R/, Rj), 

Zi,=0 m 2 '=-(h-l 2 ) m 2 "=-l' 2 k=0 

where the product of the two regular harmonics could be rewritten as a sum over a 
single harmonic, if the interest lies in the partial-wave analysis of the Fourier transform. 
It can be easily checked that (1241) reduces to the known result for the overlap of STO 
in the limit of vanishing momentum k. 

It is now interesting to compare with the related formula of Trivedi and 
Steinborn for the Fourier transform of B-function products [16J. At first glance 
the Trivedi-Steinborn result looks more compact and involves a lower number of 
summations. This is due to the favorable behavior of B-functions under the Fourier 
transform. If one is interested in STO, however, as it is often the case in quantum 
chemical or condensed matter problems, Equation (1241) provides the answer directly, 
while usage of the Trivedi-Steinborn form requires a summation over several individual 
integrals. Admittedly, for modest values of n only a small number of B-functions is 
necessary to represent a certain STO. 

There is however another point which should be important in terms of efficiency. In 
a numerical quadrature a large number of function evaluations is necessary, especially 
if one tries to achieve high precision. In the Trivedi-Steinborn form regular spherical 
harmonics appear under the pending one-dimensional integral, while the integrand in 
(122]) is simpler. Moreover, the quantity S^]^ k n212 does not depend on magnetic quantum 
numbers and can be precomputed for every l[, V 2 and stored in an array of dimension k. 



5. Implementation details 

In this paragraph we provide information on the implementation of the derived 
expressions, discuss the issue of numerical stability and give some benchmark data. 

The formulas of the last section are also valid in the special case of two STO located 
on the same center, due to the following property of the modified McDonald function 

lim x 2u LJx) = 2 v -\v - 1)! \/v > 0. (25) 

x— >0 

Nevertheless, it is computationally much more efficient to replace the STO product by 
a sum over single STO using Gaunt coefficients. In this way the known analytical result 
for the Fourier transform of individual STO given by Belkic and Taylor [35] may be 
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employed. Since the routine for the computation of real Gaunt coefficients is called 
extremely often also in the two-center case, an efficient strategy for their evaluation 
becomes very important. We follow the recent work of Pinchon and Hoggan [36], who 
devised a new index function to retrieve precomputed Gaunts for complex spherical 
harmonics. Only those coefficients that do not vanish due to selection rules are actually 
stored initially. Real Gaunt coefficients may then be obtained as outlined by Homeier 
and Steinborn [37J. 

The remaining computational bottleneck is given by the numerical integration. As 
already mentioned in the previous section, the term S^]! 1 ^ 1212 is constructed right after 
looping over l[,l' 2 as an one-dimensional temporary array. The integrals over given 
values of /x, v and a in fl23l are computed only once and then stored, since they appear 
repeatedly for different combinations of the summation variables. For the numerical 
quadrature itself, we use adaptive integration as implemented in the qag routine of the 
QUADPACK library with a (7,15) Gauss-Kronrod rule [38]. With this approach an accuracy 
of typically 14 significant figures is achieved for the basic integrals as well as the overall 
Fourier transform. 

The algorithm presented here is numerically stable for a wide range of quantum 
numbers, inter-center distances and momenta k. In situations where the ratio of decay 
constants C1/C2 is large, we however do find a significant digital erosion. For example, 
we found still 13 figure accuracy for a certain integral with a decay constant ratio of 50, 
which reduced to eleven figures at a ratio of 100 and finally three figures at a ratio of 150. 
This drawback was also observed in related earlier studies [I7J [39] and possible remedies 
were suggested by Homeier and Steinborn |40j and recently by Safouhi and Berlu [H] . In 
most real world applications the atomic numbers of elements constituting the structure 
in question usually do not differ grossly. If the interest is however in properties like 
electronic excited states or polarizabilities, additional diffuse basis functions with small 
decay constants are required. In these careful and more sophisticated evaluation 

of the basic integrals is necessary as outlined for example by Homeier and Steinborn 

In table [1] and [2] we provide some benchmark results for selected parameter 
values. The numerical error is estimated by a comparison with direct three-dimensional 
integration (Equation (JED) performed with the computer algebra package maple, that 
features arbitrary precision arithmetic. The CPU timings of the algorithm were 
performed on an Intel Pentium IV at 3.40GHz. The evaluation of a Fourier transform 
takes roughly some hundreds of /is which can be compared to the computational cost of a 
simple overlap integral on a similar machine, which was reported to be about three orders 
of magnitude lower [20]. This had to be expected, since in the latter case no numerical 
quadrature is required. In addition, Equation fl24|) shows a much higher complexity 
than the expression for the overlap. An important point for calculations in extended 
basis sets is also apparent from table [TJ The general computational cost increases 
with increasing angular momentum, but the integrals for different combinations of the 
magnetic quantum number come at little additional cost. In fact, the CPU time per 
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Table 1. Fourier transforms over products of normalized STO which share the 
following parameters: Rj = (0.3, -0.6, 0.9), R, 7 = (1.8, 0.9, 0.1), k = (0.4, -0.7, 
0.1), rii — 5, rT-2 = 4 (principal quantum number), Q\ = 3.0, £2 = 9.0. 
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Table 2. Comparison of accuracy and numerical efficiency of the algorithm presented 
in this work with the one of Trivedi and Steinborn in the implementation of Homeier 
and Steinborn [40]. Parameters for the various integrals are the same as in table [1] The 
provided number of significant digits (Digits) is the minimum of the digits for real and 
imaginary part, respectively. CPU times in ms (Time) correspond to the computation 
of (2Zi +1) x (2I2 +1) integrals and present an average over 1000 evaluations. 
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integral is decreasing with increasing I. This is a consequence of the fact that the major 
bottleneck of this scheme is the construction of the quantity S'D)! 1 ^ 1212 (122]) which is 
independent of m. 

In order to further explore the numerical efficiency of our approach, we performed 
test calculations with the FT2B code of Homeier, which implements the Trivedi- 
Steinborn formula and is described in detail in [ID] . Using Mobius-transformation-based 
quadrature rules, these authors were able to handle the highly oscillatory integrand of 
the remaining one-dimensional quadrature very efficiently. Utilizing the known formulas 
for the conversion of B-functions to STO (see e.g. [IS]), we were able to reproduce the 
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results of table [TJ with one exception^. The comparative timings given in table [2] were 
performed on the same machine and with comparable code optimization. Since the 
FT2B implementation is based on complex spherical harmonics, evaluations for different 
combinations of magnetic quantum numbers were necessary to obtain Fourier transforms 
of real STO. This additional effort was not included in the timings, since the Trivedi- 
Steinborn formula might be equally well formulated in real spherical harmonics. 

We find for the special choice of quantum numbers given in table [2], that the FT2B 
implementation is superior to our approach for individual integrals by roughly a factor of 
four. In general, one STO product may be represented by — Z) /2J + l) 2 B-function 
products, so that this result is strongly parameter dependent. In applications one is 
usually interested in the full set of integrals for different combinations of m-values and 
here our approach is numerically more efficient as table [2] shows. These computational 
savings will moreover increase with increasing angular momentum. 

Code improvements are possible for both the B-function approach as well as 
for our scheme. Homeier mentions in his dissertation [19] . that storage of some 
intermediate quantities might improve the performance for higher angular momentum. 
Our implementation might benefit from the Mobius quadrature put forward in [40J. 
While the integrand is evaluated at 36 points in the FT2B implementation, our adaptive 
integration requires 135 points for the same precision. A speed-up of a factor of four 
seems therefore achievable. 

6. Summary 

In this work a compact general purpose formula for the Fourier transform of STO 
products with arbitrary quantum numbers and geometrical parameters was derived. 
We highlighted the relation to earlier work based on B-functions and found differences 
that are relevant for the numerical efficiency. It should be stressed that the derivation 
presented here is completely independent. Moreover, the final formula can not be 
reduced to the Trivedi-Steinborn result by a mere transformation from B-functions to 
STO. Regarding numerical stability which is often an issue in STO related studies [42J, 
we achieved in general a completely satisfying accuracy apart from the known problems 
with very unsymmetric orbital products. We expect that the typical computational 
cost of several /is per integral allows for a very efficient evaluation of the notoriously 
complicated four-center electron repulsion integrals. The Fourier transform technique 
hence provides a viable alternative to existing direct methods in the field. 
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Appendix A. Some derivatives and further definitions 

The derivative of the modified McDonald function k v (x) has the following simple form 

= -xk v -i{x). (A.l) 

dx 

-.nil 



In order to evaluate the quantity h%££ in ([To]) an expression for the repeated action of 

ld_ 

k v ( RjQ + k 2 u(l-u)) = R 2l u%^ ( RJQ + k 2 u(l - u) 



the operator — -Jp: on k is required. Straightforward differentiation leads to 
/ 1 d x ; 



Ci dCiJ 
1 8 N ; 



K ( i2 VC^ + Wu{l - u) ) = R 2l (l - u)%-i ( i? JQ + fc 2 «(l - w) ) . (A.2) 



The action of the operator — is more involved but can be reduced to (1A.2I) . 

/ fl\"_ n!(-l)" A (-2C 2 ) j (I d\ 

\ d() " (2C)« . = ^{2i-n)\{n-i)\ \ ( d( ) ' 1 ' ) 

which gives rise to the definitions of the coefficients cf(Q m dSSJ). 

We now prove (1A. 3|) by using induction. The induction basis for n — 1 is trivial. 
We further have (5^ denoting the Kronecker delta) 

(_d\(_d_Y v n!(-l)"+^2— w ^_IAV 

I 9CA ^(2z-n-l)!(n-.)! C V 1 «W»J^ ( g ( ) 

where we used the induction hypothesis and changed the summation index to i' — i + 1 
in the second sum. Separating the term for the lower limit % = L 2 ^] m the first sum 
and the upper limit i' — n + 1 in the second sum, we arrive at 

(_d_Y +1 n!(-irW + i 2 L^J-» 2 |^,_ B _w ^M^V^ 

n 



(n + l)!(-l)"+W2— 1 _ t /ljT 
^ (2i-n- l)!(n+l-i)! \ C^C, 



i= L 2±ij+i 

.n+1 ' 1 ^ X 



For even n, we have L^pJ = n /2, and the first term in (1A.5|) vanishes. Since in this 
case L^J + 1 = L^J , ^ follows: 

(_d_Y +1 = "f 1 (» + i)!(-i)" +1+i 2-- 1 c2M ( idV 

{ d() fa {2 i-n-l)\(n + l-i)^ \ (d() ' 1 ] 

% ~ L 2 J 

For odd n, we have [^J = (n + l)/2, as well as [^J + 1 = L^J + L If we now 
extend the second sum in (1A.5|) to the lower limit i = [^f^J the compensating term 
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cancels exactly the first term in (IA.50 . Also in this case we therefore arrive at (1A.6[) . 
that is the hypothesis for n + 1, which was to be demonstrated. 

Derivatives with respect to the inter-center distance R are likewise readily obtained: 

[km) R ~ 2V " K (^ 2 + fc Ml-«)) = (-l) l R- 2il/+l) k + l (Ry/£ + k*u(l-u?) , (A.7) 
where we have used the following recursion: 

k v +i(x) = 2vk v (x) + xk„-\{x). (A. 8) 

A combination of (I A. 21) to (1A.7|) leads to second line in (122]) . 



Finally, we define the coefficients c^™* 21712 that appear in ( 12T1) . A product of two 
regular harmonics with same argument can be linearized as follows: 

z-(R)C(R) ^ECCW^ 2 "" ( A -9) 

In terms of the Gaunt-like coefficients for the real unnormalized spherical harmonics of 

[hm 1 \l 2 m 2 \l z m z ] = J Y hmi {n)Y hm2 {n)Y l3m3 (Q)dQ, (A.IO) 
these coefficients read 

Please note that notation (lA.lOj) differs from the one usually employed for Gaunt 
coefficients [43J. The linearization formula f]A.9|) may be considerably simplified by 
taking advantage of the selection rules for the Gaunt-like coefficients ( IA.10I) . which were 
discussed by Homeier and Steinborn |37j : 

^7(R)C(R) =EEteW)^ 

Z=0 m 

m G {mi + TTi2, m>i — wii-, —m-i + —mi — m-i\ (A. 12) 

Appendix B. Leibniz theorem for regular harmonics 

Here we prove Eq. ([171) of section @] 

^ m (V)(/s) = E £ £ d l r m , m „ (zfAV) f) (zf (V) g) . (B.l) 

l'=0 m'=-(l-l') m"=-V 

The regular harmonic ^["(V) is given in cartesian form as 

l l-i ,Q\W ,Q\U) ,Q\<!-i-3) 



where (d/dx)^ denotes the n-th partial derivative with respect to x and the coefficients 



i=0 j=0 



C^J 1 are known constants (see e.g. [29], Eqs. 3 and 4). Applying the Leibniz theorem 
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for the differentiation of products we have 

*T(v) (/ 9 ) = E E eg" E E E J (•,) ( fl ~ i ~ j 

i=0 j=0 i'=0 j'=0 k'=0 

(*-*') / a \ (jWO 



X 



(B-3) 



<9x 



9_Y° (jlY' (jl 



(k>) 



dx 



dy I \dz 



9 



Using now the completeness of the regular harmonics we may expand the product x l y 3 z k 
into harmonics of angular momentum I = i + j + k 



x i y j z k 



E B l m zf{v) ; l = i+j + k 

m=—l 

B l m = J x i y>z k r- 2l zY l (T) dQ. 
Inserting this expansion into (IB. 31) we find 



(B.4) 



x 



i'+j'+k' 



B i'+r+k'm" f V \ 
/ . n m" z i'+j'+k'\ w 1 9 



^ ^(v)(/ 9 )=^:I:c'^■i:^:'i: ^ (;,)(")( , 1 3 

i=0j=0 i'=0 j'=0 k'=0 V / V/ \ h 

1 E " B l J'-^ y zT: v ^ kl {v)f 

_ m '=-(l-i'-j'-k') J [«■»"=-(*'+.?'+*') 

which can be simplified after changing the summation order according to 

a' b 1 c' a'+b'+c' min(a,6') min(a— 6,c') 

EEEf(4c)= EE E F(a-6-c,6 )C ), (B.6) 

a=0 6=0c=0 a=0 fe=0 c=max(0,a-b-a') 

for arbitrary F. With the help of the coefficients A l [ m 

l i-i min(Z'j) mm(l'-j',l-i-j) / • \ / • \ //_.•_ A 

i=0j=0 j'=0 fc'=max(0,Z'-j'-i') V J ^ / V / \ h I 

with I' — i' + j' + k' , we finally arrive at 

i i-i' v 



(B.5) 



l (v)(/^) = E< E E K~^-V(v)/ 

l'=0 m'=-(l-V) m"=-V 



Bl„zf{V)g ,(B.8) 



which is equivalent to flUD if we set rf]^, m „ = A^B^B^, 1 . 



Appendix C. Proof of equation ( 119ft 

An old theorem given by Hobson [51] (see also [27]) states that if R = (X 2 + Y 2 + Z 2 )^ 
and H(X, Y, Z) is a homogeneous polynomial of degree I in the X, Y, Z variables, then: 

( 8 8 8 \ 1 2 l - 2k r i ( 8 \ l ~ k 

Taking into account that -J^ = — = from R = Rj — R/, we may apply (IC.ll) 
to the product of the regular harmonics ^ X (V) and z^ 2 (V), which is a homogeneous 
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polynomial of degree li + l 2 . Hence, 



h+h 2h+h-2k 



Z, 



r(V/)C(Vj)/(i2) = (-i) 11 E 

fc=0 



fc! 



-i)' 1 E 



k=0 



k\ 



v 2fe ^(R)C(R-) 



/(*) 



/(i2),(C2) 



which is the same as (1191) apart from the upper limit in the sum over k. In Ref. poj it 
was shown that 



V 2k z™(R)R 



21 



2 2k l\F( P +l+3/2) , m /p^p2l-2fe . h<] 
(l-k)\r(p+l-k+3/2) Z p \ JX )- n ' ■ h 1 



(C.3) 



: k> I 

Combining fl02l) with (jXl2l> and ({C~3]) we arrive at p5|), (|20j) and of the main 
paper. 
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